n <- 50 # number of trajectories per noise level
noises <- seq(0.2, 3, 0.2)
# Phase shifted
ps <- multi_sims(type = "ps", noises = noises, n = n)
ps
## Time variable value noise
## 1: 0 V1 0.0000000 0
## 2: 1 V1 0.8414710 0
## 3: 2 V1 0.9092974 0
## 4: 3 V1 0.1411200 0
## 5: 4 V1 -0.7568025 0
## ---
## 79996: 95 V50 -0.9858802 3
## 79997: 96 V50 -0.6735795 3
## 79998: 97 V50 0.2580071 3
## 79999: 98 V50 0.9523832 3
## 80000: 99 V50 0.7711425 3
# Phase shifted with trend
pst <- multi_sims(type = "pst", noises = noises, n = n, slope = 0.1)
# Amplitude noise
na <- multi_sims(type = "na", noises = noises, n = n)
plot_sim(ps)
plot_sim(pst)
plot_sim(na)
For each condition, this computes the distance between the average trajectory and each individual trajectory.
cond <- "noise"
ti <- "Time"
mea <- "value"
lab <- "variable"
ps_mean <- dist_mean(data = ps, condition = cond, tcol = ti, measure = mea, label = lab, return.mean = F)
ps_mean
## noise variable euclid_to_mean
## 1: 0 V1 0.000000
## 2: 0 V2 0.000000
## 3: 0 V3 0.000000
## 4: 0 V4 0.000000
## 5: 0 V5 0.000000
## ---
## 796: 3 V46 6.928511
## 797: 3 V47 6.365261
## 798: 3 V48 7.451127
## 799: 3 V49 6.951201
## 800: 3 V50 7.809669
pst_mean <- dist_mean(data = pst, condition = cond, tcol = ti, measure = mea, label = lab, return.mean = F)
na_mean <- dist_mean(data = na, condition = cond, tcol = ti, measure = mea, label = lab, return.mean = F)
For each condition, take each pair of trajectories and compute the correlations between them (Pearson, Spearman and Kendall). In addition it clips the trajectories by comparison with their rolling means: whenever the trajectory is above its rolling mean, set it to 1 otherwise to 0. Finally for each pair of trajectories, the “overlap” metric, represents how often the trajectories of the pair have same value (i.e. 0.5 if complete random, 1 if they fully overlap).
This takes quite a while (partially) because it is implemented with for loops, but is still performed in reasonable time (i.e. a few minutes)
ps_pw <- all_pairwise_stats(data = ps, condition = cond, label = lab, measure = mea, k_roll_mean = 5)
ps_pw
## noise Label1 Label2 Overlap Pearson Spearman Kendall
## 1: 0 V1 V2 1.00 1.0000000 1.0000000 1.0000000
## 2: 0 V1 V3 1.00 1.0000000 1.0000000 1.0000000
## 3: 0 V1 V4 1.00 1.0000000 1.0000000 1.0000000
## 4: 0 V1 V5 1.00 1.0000000 1.0000000 1.0000000
## 5: 0 V1 V6 1.00 1.0000000 1.0000000 1.0000000
## ---
## 19596: 3 V47 V49 0.60 0.3140454 0.3281368 0.2222222
## 19597: 3 V47 V50 0.04 -0.9959261 -0.9967597 -0.9539394
## 19598: 3 V48 V49 0.08 -0.9488995 -0.9511191 -0.8129293
## 19599: 3 V48 V50 0.72 0.6674754 0.6372277 0.4553535
## 19600: 3 V49 V50 0.36 -0.3983772 -0.3933633 -0.2682828
pst_pw <- all_pairwise_stats(data = pst, condition = cond, label = lab, measure = mea, k_roll_mean = 5)
na_pw <- all_pairwise_stats(data = na, condition = cond, label = lab, measure = mea, k_roll_mean = 5)
cor(ps_pw[, 4:7])
## Overlap Pearson Spearman Kendall
## Overlap 1.0000000 0.9888158 0.9905337 0.9983643
## Pearson 0.9888158 1.0000000 0.9997799 0.9903528
## Spearman 0.9905337 0.9997799 1.0000000 0.9920946
## Kendall 0.9983643 0.9903528 0.9920946 1.0000000
cor(pst_pw[, 4:7])
## Overlap Pearson Spearman Kendall
## Overlap 1.0000000 0.9893778 0.9896651 0.9855989
## Pearson 0.9893778 1.0000000 0.9998151 0.9731459
## Spearman 0.9896651 0.9998151 1.0000000 0.9743867
## Kendall 0.9855989 0.9731459 0.9743867 1.0000000
cor(na_pw[, 4:7])
## Overlap Pearson Spearman Kendall
## Overlap 1.0000000 0.9527990 0.9539770 0.9668563
## Pearson 0.9527990 1.0000000 0.9963428 0.9805116
## Spearman 0.9539770 0.9963428 1.0000000 0.9806409
## Kendall 0.9668563 0.9805116 0.9806409 1.0000000
For each condition, count the proportion of trajectories that are lying more than x standard deviation from the mean (conversion to z-score first?).
Get a p-value: Assume that they are normally distributed, check that the distributions lies above 0 for correlations and above 0.5 for clipping (t-test).
ps_pw[noise !=0 , .(p.val.pears = t.test(Pearson)$p.value,
p.val.spear = t.test(Spearman)$p.value,
p.val.kenda = t.test(Kendall)$p.value,
p.val.ovrlp = t.test(Overlap, mu = 0.5)$p.value), by = noise]
## noise p.val.pears p.val.spear p.val.kenda p.val.ovrlp
## 1: 0.2 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## 2: 0.4 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## 3: 0.6 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## 4: 0.8 1.690953e-215 3.200387e-218 4.495175e-211 1.056645e-210
## 5: 1.0 1.066447e-101 7.696656e-101 3.318720e-97 4.142795e-94
## 6: 1.2 2.225065e-36 1.898027e-36 9.501207e-36 6.581236e-37
## 7: 1.4 5.429424e-12 3.994100e-12 4.904341e-12 2.852469e-13
## 8: 1.6 1.797949e-01 1.893047e-01 1.820744e-01 2.554369e-01
## 9: 1.8 2.359007e-01 2.308119e-01 2.753541e-01 2.952123e-01
## 10: 2.0 5.539425e-01 5.743377e-01 5.924177e-01 6.287601e-01
## 11: 2.2 3.331276e-01 3.178592e-01 2.714869e-01 2.239812e-01
## 12: 2.4 6.928466e-01 6.932043e-01 6.288130e-01 5.919235e-01
## 13: 2.6 4.498626e-01 4.578279e-01 5.315049e-01 5.422259e-01
## 14: 2.8 6.977451e-01 7.108618e-01 7.185016e-01 5.959808e-01
## 15: 3.0 6.402742e-01 6.340054e-01 5.474272e-01 5.356297e-01
pst_pw[noise !=0 , .(p.val.pears = t.test(Pearson)$p.value,
p.val.spear = t.test(Spearman)$p.value,
p.val.kenda = t.test(Kendall)$p.value,
p.val.ovrlp = t.test(Overlap, mu = 0.5)$p.value), by = noise]
## noise p.val.pears p.val.spear p.val.kenda p.val.ovrlp
## 1: 0.2 0 0 0 0.000000e+00
## 2: 0.4 0 0 0 0.000000e+00
## 3: 0.6 0 0 0 0.000000e+00
## 4: 0.8 0 0 0 9.563845e-188
## 5: 1.0 0 0 0 1.317767e-108
## 6: 1.2 0 0 0 1.151031e-37
## 7: 1.4 0 0 0 7.837262e-14
## 8: 1.6 0 0 0 4.045208e-01
## 9: 1.8 0 0 0 1.321593e-02
## 10: 2.0 0 0 0 4.791700e-02
## 11: 2.2 0 0 0 7.028294e-01
## 12: 2.4 0 0 0 5.277985e-03
## 13: 2.6 0 0 0 3.453810e-01
## 14: 2.8 0 0 0 8.824005e-01
## 15: 3.0 0 0 0 9.641060e-01
na_pw[noise !=0 , .(p.val.pears = t.test(Pearson)$p.value,
p.val.spear = t.test(Spearman)$p.value,
p.val.kenda = t.test(Kendall)$p.value,
p.val.ovrlp = t.test(Overlap, mu = 0.5)$p.value), by = noise]
## noise p.val.pears p.val.spear p.val.kenda p.val.ovrlp
## 1: 0.2 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## 2: 0.4 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## 3: 0.6 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## 4: 0.8 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## 5: 1.0 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## 6: 1.2 0.000000e+00 0.000000e+00 0.000000e+00 6.583834e-294
## 7: 1.4 0.000000e+00 0.000000e+00 0.000000e+00 1.110489e-183
## 8: 1.6 0.000000e+00 0.000000e+00 0.000000e+00 2.473046e-139
## 9: 1.8 5.115229e-298 1.423070e-275 1.753057e-274 1.609504e-70
## 10: 2.0 7.239534e-229 7.302942e-219 1.839038e-217 8.695974e-72
## 11: 2.2 2.623503e-154 5.236024e-153 1.867011e-151 4.982420e-36
## 12: 2.4 1.425185e-88 5.410524e-82 9.524339e-81 4.832441e-14
## 13: 2.6 9.667281e-95 6.927798e-90 4.912993e-90 5.627369e-29
## 14: 2.8 2.233522e-102 7.159980e-97 9.882664e-97 1.278583e-23
## 15: 3.0 4.504603e-77 3.995172e-79 1.910595e-78 7.032621e-24